Hyperacuity system and methods for real time and analog detection and kinematic state tracking

ABSTRACT

Certain embodiments of the methods and systems disclosed herein determine a location of a tracked object with respect to a coordinate system of a sensor array by using analog signals from sensors having overlapping nonlinear responses. Hyperacuity and real time tracking are achieved by either digital or analog processing of the sensor signals. Multiple sensor arrays can be configured in a plane, on a hemisphere or other complex surface to act as a single sensor or to provide a wide field of view and zooming capabilities of the sensor array. Other embodiments use the processing methods to adjust to contrast reversals between an image and the background.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application is a Divisional of U.S. patent application Ser. No.14/935,352 filed on Nov. 6, 2015, which is incorporated herein byreference. U.S. patent application Ser. No. 14/935,352 claims thebenefit of U.S. Provisional Patent Application No. 62/077,095, filedNov. 7, 2014, and U.S. Provisional Patent Application No. 62/079,688,filed Nov. 14, 2014, each of which is incorporated herein by reference.

BACKGROUND OF THE INVENTION

Optical sensing, detecting, and tracking devices may comprise a lenssystem to project light rays from a segment of an environment onto aplanar or other surface. One such device is shown schematically inFIG. 1. The surface may comprise an array of photo sensors that createelectronic signals from the projected light. Examples of photo sensorsinclude charge coupled devices (CCD) having a planar rectangular arrayof distinct light receptors that produce a respective pixel of theprojected light, such as are used in digital cameras. Typically eachsuch receptor in a CCD generates a voltage corresponding to the lightamplitude level integrated over the exposure time. Also, the values ofthe pixels in a CCD are typically taken at specific time instants,rather than having a continuous time, or real-time, output. That is, thereceived light is digitally sampled in time by the CCD rather thanproducing a real-time analog output. Such devices typically includeelectronics that may amplify the sensor signals from the optical sensorsand applying signal processing to the sensor signals to obtain desiredoutput results, such as position, motion and acceleration of an object.Other types of photo sensors may be chosen, depending on theapplication.

Depending on the application, other detecting and tracking devices maybe based on non-visible electromagnetic radiation. Still other detectingand tracking devices may be based on non-electromagnetic phenomena.

However, many current devices using multiple sensors typically areconfigured so that the sensors only react or respond to inputs local tothemselves. Inputs impinging on other nearby sensors are considered tobe problems to be avoided or designed around. This limits the overallresolution of, for example, location detection.

BRIEF SUMMARY OF THE INVENTION

Systems and methods are disclosed for object detection and trackingusing multiple sensors, such as photo sensors or acoustic wave sensors,that are arranged to have overlapping output response profiles. Thesensors' outputs can then be differenced or combined to produce alocation accuracy at least, and often better, than what could beachieved by the same sensors without overlapping profiles. Anarrangement of sensors having such overlapping output response profilesis termed an ommatidium (plural: ommatidia).

Hyperacuity is the ability of a detection or tracking device or systemto discriminate a minimal angular displacement of an object moving inthe system's field of view that is smaller than the sensor dimension orspacing. The hyperacuity of such systems can be limited only by thesystem electronics. In preferred embodiments, analog electronics areused to achieve faster results by avoiding limitations due to samplingrate; nevertheless, the methods and systems can still be implemented indigital electronics.

Multiple ommatidia can be organized into a single system to achievevarious advantages, such as increased resolution and/or wider field ofview. Also, a multiple ommatidia can be organized hierarchically toachieve zooming capability for detection and tracking.

Multiple ommatidia can be configured on surface to provide a wide fieldof view (FoV) for detection and tracking. The multiple ommatidia can usecombined methods of detection and tracking, and may use IR, visible, orpolarization detection.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 shows a structure of a tracking system comprising an ommatidium,according to an embodiment.

FIGS. 2A, 2B, 2C, and 2D shows configurations of sensors arranged toform a single ommatidium to achieve hyperacuity, according toembodiments.

FIG. 3 shows an embodiment of a sensor configured with a lens, accordingto embodiments.

FIG. 4 shows overlapping Gaussian response profiles of sensor signalswithin an ommatidium, according to an embodiment.

FIGS. 5A and 5B respectively show a functional diagram for a circuit toimplement a differencing method and Tables of coefficients to use withthe circuit components, according to embodiments.

FIGS. 6A and 6B respectively show the geometric parameters underlying anorthogonal decomposition method and a functional diagram for a circuitto implement an orthogonal decomposition method, according toembodiments.

FIGS. 7A, 7B and 7C shows the geometric parameters used in an analysisof an orthogonal decomposition method, according to an embodiment.

FIG. 8 shows graphs related to a sensor output during traversal of anoptical sensor by an object image, according to embodiments.

FIG. 9 shows the geometry used for analyzing a trajectory of a movingobject's image, according to embodiments.

FIG. 10 shows a functional block diagram of operations used to track amoving object image, according to embodiments.

FIG. 11 shows a component level block diagram of operations used totrack a moving object image, according to embodiments.

FIGS. 12A and 12B show the geometry used for determining a trajectory ofa moving object image, according to embodiments.

FIG. 13 shows an component level block diagram of operations used totrack a moving object image, according to embodiments.

FIG. 14 shows a process of dynamic contrast inversion, according toembodiments.

FIG. 15 shows a component level block diagram of operations used toimplement dynamic contrast inversion, according to embodiments.

FIGS. 16A, 16B and 16C show configurations of sensors and structures formulti-ommatidia systems and methods, according to embodiments.

FIGS. 17A, 17B and 17C show the geometry of a multi-ommatidia system anda corresponding table of parameters, according to embodiments.

FIG. 18 shows the geometry of a multi-ommatidia system, according toembodiments.

FIG. 19 shows the geometry for tracking an object image betweenommatidia of a multi-ommatidia system, according to embodiments.

FIG. 20 shows a zooming method of a multi-ommatidia system, according toembodiments.

FIG. 21 shows a zooming method of a multi-ommatidia system, according toembodiments.

FIG. 22 is a flow chart of a method of tracking an object image in amulti-ommatidia system, according to embodiments.

FIG. 23 shows a multi-ommatidia system having adaptive grouping ofsensors into ommatidia, according to embodiments.

FIG. 24 shows a functional block diagram for adaptive grouping ofsensors into ommatidia within a multi-ommatidia system, according toembodiments.

FIG. 25 is a block diagram for using orthogonal differencing methods ina multi-ommatidia system, according to embodiments.

FIG. 26 is a block diagram for combined methods in a multi-ommatidiasystem, according to embodiments.

FIG. 27 is a component block diagram for combined methods in amulti-ommatidia system, according to embodiments.

DETAILED DESCRIPTION OF THE INVENTION I. Introduction

Systems and methods are disclosed in which multiple sensors are arrangedto form a combined sensing unit, called an ommatidium. The sensors maybe of any of a variety of physical phenomena, such as visible, infraredor other electromagnetic waves, acoustic waves, or light polarization.

In exemplary embodiments the sensors in an ommatidium have overlappingoutput response profiles. This overlapping response is typicallynon-linear, and is in fact used to achieve a higher location resolution(hyperacuity) than would a system using similar sensors without suchoverlap. Further, exemplary embodiments will make use of analogcomponents, such as fast logarithmic amplifiers, to achieve fastdetection and tracking that is not limited by a sampling rate, as in adigital system. The methods and systems can nevertheless be implementedwith sampling and digital processing.

For simplicity of exposition, the exemplary embodiments disclosed hereinwill be of methods for, and systems of, photo sensors but it will beclear to one of skill in the art that the methods and systems are notrestricted to photo sensor systems. The methods and systems describedherein for optical (visual light) based sensing and tracking using photosensors will be seen as readily adaptable to other frequency ranges(such as infrared, ultraviolet or x-ray) of electromagnetic radiationfor corresponding sensors and imaging devices. Further, the methods andsystems in fact will be seen as readily adaptable to location andtracking devices based on other physical phenomena, such as for sonardevices that sense acoustic waves.

The following sections describe different embodiments of the inventions.Part II discloses systems and methods based on a single ommatidium arrayof sensors. Methods of achieving hyperacuity include differencingmethods that may use logarithmic amplification of sensor signals andanalog electronics for real-time tracking. Other methods disclosedinclude orthogonal decomposition methods that use sinusoidal weightingof the sensor signals to form a pair of resultant signals which are usedto achieve hyperacuity tracking. Part III discloses systems and methodsbased on multiple ommatidia arrays. Such systems may use a single lensor multiple lenses. Part IV discloses embodiments directed to HYPRIS'systems, which may use multiple ommatidia systems to achieve wide fieldof view, and other tracking and imaging results.

II. Single Ommatidia Systems and Methods

The methods and systems disclosed in this section combine multiplesensors to function as a single location or position detection (hereinjust “detection”) and location tracking (herein just “tracking”) device.Part A describes exemplary physical system structures that may be usedto implement the systems and methods. Part B describes methods based ondifferencing, and Part C describes alternate methods based on orthogonaldecomposition of the sensor signals. Part D describes how the methodsand systems can implement trajectory tracking. Part E describes methodsto compensate when a sensor's output signal response to an image featureof a detected/tracked object in fact is darker than the output responseto the imaged background. Part F describes ways to combine the variousmethods.

A. Single Ommatidium Systems

FIG. 1 shows an exemplary detection and tracking system 100 that canmake use of the methods and systems disclosed herein. In this example, asupporting structure 120 encloses a lens 110, with a focal length f,140, that focuses light onto a sensor field 130. The sensor field cancomprise a plurality of sensor arrangements, as now described.

FIGS. 2A, 2B, 2C, and 2D show exemplary configurations of sensors:respectively 210, 220, 230 and 240. For the case that the sensors arephoto sensors, these can be arranged on the sensor field 130. FIG. 2Ashows a first embodiment of an ommatidium configuration, in which sevensensors are arranged around a central sensor. Depending on the shape ofthe individual sensors, a similar configuration is shown in FIG. 2B forcircular sensors.

In FIG. 2B the sensors R₀, . . . R₆ are arranged geometrically to forman ommatidium disposed in a hexagonal pattern around the central sensorR₀. Each sensor has a unique and well-defined location relative to thecentral sensor which is also the central or reference point of theommatidium. When composed of photo sensors, a beam of light incident onthis array of sensors will form an image spot on a given sensor,indicated at point (x, y) on R₃. Depending the output response profileof the sensors, as shown in FIG. 3, one or more of the lenses may usethe configuration 300 in which there is a surface mounted lens 310 onthe sensor. The combination of lens 110 and lens 310 may implement aresponse profile of the sensor that is Gaussian and/or overlaps with theresponse profile of neighboring sensors. An example of overlappingresponse profiles for the configuration of FIG. 2B is shown in FIG. 4and described further below.

The described sensor response profile is to be a monotonicallydecreasing function from each sensor's center; these decaying functionsinclude but are not limited to Gaussian, Airy, cosine power, sin(x)/x,sinc(x), (1−x²)^(n) polynomials, etc. There are several ways to generatethe nonlinear response profile, including nonlinear doping of thedetector, optical termination, waveguide transmission, and intentionalimage feature blurring to achieve the desired response profile. Eachsensor has an overlapping field of view with its adjacent neighbors.

FIGS. 2C and 2D show ommatidia with alternate arrangements 230 and 240of sensors. The methods described below can be also be adapted to theseconfigurations, as will be clear to one of skill in the art afterunderstanding the following descriptions.

B. Differencing Systems and Methods

This part describes differencing methods by which the exemplaryarrangements of multiple individual sensors into ommatidia previouslydescribed can implement hyperacuity detection and tracking.

Embodiments of the differencing method will be disclosed by way ofexample using the seven sensor hexagonal ommatidium configuration ofFIG. 2B. From the derivation given below it will be clear to one ofskill in the art how the methods can be adapted for ommatidia havingother configurations.

For this example case, each sensor is assumed to have a Gaussian-type,axisymmetric response profile, as exemplified in FIG. 4, that overlapswith adjacent sensors. An image spot is shown located at position (x, y)in the array's coordinate system in FIG. 2B. An image feature of theobject produces a light spot that is sensed by several neighboringsensors due to their overlapping fields of view. The location of theimage spot's center as observed from the center of the photo sensors R₀,R₂ and R₃ in FIG. 2B can be expressed by the set of equations

(x-x ₀)²+(y-y ₀)² =r ₀ ²

(x-x ₂)²+(y-y ₂)² =r ₂ ²

(x-x ₃)²+(y-y ₃)² =r ₃ ²  (II.B.1)

The solution of this set of equations can be interpreted as theintersection of three circles at the point (x, y), the circles havingradii r₀, r₂ and r₃ centered at (x₀, y₀), (x₂, y₃), and (x₃, y₃),respectively. Taking the difference between the second and the firstequation in (II.B.1), and then taking the difference between the thirdand the first in (II.B.1) yields

(x-x ₂)²+(y-y ₂)²−(x-x ₀)²−(y-y ₀)² =r ₂ ²-r ₀ ²

(x-x ₃)²+(y-y ₃)²−(x-x ₀)²−(y-y ₀)² =r ₃ ²-r ₀ ²  (II.B.2)

respectively. Letting the point (x₀, y₀)=(0, 0) be the origin, whichcoincides with the geometric center of the array, the set (II.B.2) thenbecomes

(−2x ₂)x+(−2y ₂)y=r ₂ ²-r ₀ ² −x ₂ ²-y ₂ ²

(−2x ₃)x+(−2y ₃)y=r ₃ ²-r ₀ ² −x ₃ ² y ₃ ²  (II.B.3)

The right hand sides of (II.B.3) are constant values:

A=r ₂ ²-r ₀ ² −x ₂ ²-y ₂ ²

B=r ₃ ²-r ₀ ² −x ₃ ²-y ₃ ²  (II.B.4)

Using Cramer's rule one finds from (II.B.3) and (II.B.4) the unknowns xand y as

$\begin{matrix}{{x = {\frac{1}{2}\left( \frac{{By}_{2} - {Ay}_{3}}{{x_{2}y_{3}} - {x_{3}y_{2}}} \right)}},{y = {\frac{1}{2}\left( \frac{{Ax}_{3} - {Bx}_{2}}{{x_{2}y_{3}} - {x_{3}y_{2}}} \right)}}} & \left( {{{II}.B}{.5}} \right)\end{matrix}$

Once the x and y coordinates of the image-spot have been determined, theangle φ that the spot subtends with respect to the array's referenceaxis and the distance ρ to the origin and center of the array are givenrespectively by

$\begin{matrix}\begin{matrix}{{\phi = {\arctan \left( \frac{y}{x} \right)}},\mspace{14mu} {\rho = \left( {x^{2} + y^{2}} \right)^{1/2}}} & \;\end{matrix} & \left( {{{II}.B}{.6}} \right)\end{matrix}$

The term r_(i) in these equations, i.e. the distances from the imagespot center to the center of the photoreceptors, is determined byanalyzing the response profile of the ommatidium's sensors. The Gaussiansensor response profile to focused light stimuli has the form:

I _(i) =Îe ^(−βr) ²   (II.B.7)

Equation (II.B.7) represents a Gaussian-like radial-dependence response.In this equation I_(i)=I(r_(i)) is the signal intensity or response ofthe i^(th) sensor to an image spot illuminating the sensor at a distancer_(i) from the sensor's axis of symmetry, β is a constant thatrepresents the rate of decay of the response with the distance from thecenter of the sensor. One may relate β to the value σ of the Gaussiansused in statistics, the conversion is β=1/(2σ²), and Î is the peak valueof the sensor response for a given image spot illuminating the center ofthe sensor. A cylindrical coordinate system is used here for the sensorresponse function.

The square of the radial distance at which the image-spot isilluminating the sensor, relative to the sensor's center, can beobtained from (II.B.7) as

$\begin{matrix}{r_{i}^{2} = {\frac{1}{\beta}{\ln\left( \frac{\hat{I}}{I_{i}} \right)}}} & \left( {{{II}.B}{.8}} \right)\end{matrix}$

Equation (II.B.8) requires knowing Î to solve for r_(i) since only I_(i)the sensor signal response due to the image spot illuminating at adistance r_(i) is observed or measured. It will be shown below that theknowledge of this value is not needed, as relative measures of thesensor responses will only be required.

Referring back to FIG. 2B and using equation (II.B.8), one has

${r_{2}^{2} - r_{0}^{2}} = {{{\frac{1}{\beta}{\ln\left( \frac{\hat{I}}{I_{2}} \right)}} - {\frac{1}{\beta}{\ln\left( \frac{\hat{I}}{I_{0}} \right)}}} = {\frac{1}{\beta}\left\lbrack {{\ln\left( \frac{\hat{I}}{I_{2}} \right)} - {\ln\left( \frac{\hat{I}}{I_{0}} \right)}} \right\rbrack}}$

which becomes

$\begin{matrix}{{r_{2}^{2} - r_{0}^{2}} = {\frac{1}{\beta}{\ln \left( \frac{I_{0}}{I_{2}} \right)}}} & \left( {{{II}.B}{.9}} \right)\end{matrix}$

As shown in (II.B.9), Î was canceled-out by a logarithmic property.Similarly,

$\begin{matrix}{{r_{3}^{2} - r_{0}^{2}} = {\frac{1}{\beta}{\ln\left( \frac{I_{0}}{I_{3}} \right)}}} & \left( {{{II}.B}{.10}} \right)\end{matrix}$

Substituting (II.B.9) and (II.B.10) into the right-hand side of (4)produces

$\begin{matrix}{{A = {{{\frac{1}{\beta}{\ln \left( \frac{I_{0}}{I_{2}} \right)}} - \left( {x_{2}^{2} + y_{2}^{2}} \right)} = {{\frac{1}{\beta}{\ln \left( \frac{I_{0}}{I_{2}} \right)}} + d^{2}}}}{B = {{{\frac{1}{\beta}{\ln \left( \frac{I_{0}}{I_{3}} \right)}} - \left( {x_{3}^{2} + y_{3}^{2}} \right)} = {{\frac{1}{\beta}{\ln \left( \frac{I_{0}}{I_{3}} \right)}} + d^{2}}}}} & \left( {{{II}.B}{.11}} \right)\end{matrix}$

In (II.B.11), d is the constant distance from the center of the sensorarray to each photoreceptor. Substituting A and B from (II.B.11) into(5) give the image spot coordinates (x, y) as

$\begin{matrix}{x = \frac{{\frac{y_{2}}{\beta}{\ln \left( \frac{I_{0}}{I_{3}} \right)}} - {\frac{y_{3}}{\beta}{\ln \left( \frac{I_{0}}{I_{2}} \right)}} + {d^{2}\left( {y_{3} - y_{2}} \right)}}{2\left( {{x_{2}y_{3}} - {x_{3}y_{2}}} \right)}} & \left( {{{II}.B}{.12}} \right) \\{y = \frac{{\frac{x_{3}}{\beta}{\ln \left( \frac{I_{0}}{I_{2}} \right)}} - {\frac{x_{2}}{\beta}{\ln \left( \frac{I_{0}}{I_{3}} \right)}} - {d^{2}\left( {x_{3} - x_{2}} \right)}}{2\left( {{x_{2}y_{3}} - {x_{3}y_{2}}} \right)}} & \left( {{{II}.B}{.13}} \right)\end{matrix}$

For the general case, one would have

$\begin{matrix}{x = \frac{{\frac{y_{i}}{\beta}{\ln \left( \frac{I_{0}}{I_{j}} \right)}} - {\frac{y_{j}}{\beta}{\ln \left( \frac{I_{0}}{I_{i}} \right)}} + {d^{2}\left( {y_{j} - y_{i}} \right)}}{2\left( {{x_{i}y_{j}} - {x_{j}y_{i}}} \right)}} & \left( {{{II}.B}{.14}} \right) \\{y = \frac{{\frac{x_{j}}{\beta}{\ln \left( \frac{I_{0}}{I_{i}} \right)}} - {\frac{x_{i}}{\beta}{\ln \left( \frac{I_{0}}{I_{j}} \right)}} - {d^{2}\left( {x_{j} - x_{i}} \right)}}{2\left( {{x_{i}y_{j}} - {x_{j}y_{i}}} \right)}} & \left( {{{II}.B}{.15}} \right)\end{matrix}$

Applying some algebra, (II.B.14) and (II.B.15) can be written as

$\begin{matrix}{x = \frac{{\left( {y_{j} - y_{i}} \right)\left\lbrack {d^{2} - {\frac{1}{\beta}{\ln \left( I_{0} \right)}}} \right\rbrack} + {\frac{1}{\beta}\left\lfloor {{y_{j}{\ln \left( I_{i} \right)}} - {y_{i}{\ln \left( I_{j} \right)}}} \right\rfloor}}{2\left( {{x_{i}y_{j}} - {x_{j}y_{i}}} \right)}} & \left( {{{II}.B}{.16}} \right) \\{y = \frac{{\left( {x_{j} - x_{i}} \right)\left\lbrack {{\frac{1}{\beta}{\ln \left( I_{0} \right)}} - d^{2}} \right\rbrack} - {\frac{1}{\beta}\left\lbrack {{x_{j}{\ln \left( I_{i} \right)}} - {x_{i}{\ln \left( I_{j} \right)}}} \right\rbrack}}{2\left( {{x_{i}y_{j}} - {x_{j}y_{i}}} \right)}} & \left( {{{II}.B}{.17}} \right)\end{matrix}$

Since the only variables in (II.B.16) and (II.B.17) are the signals I₀,I_(i), and I_(j), these equations can be written, respectively, as

x=η _(i) ln(I _(i))+η_(j) ln(I _(j))+η₀ ln(I ₀)+η_(d)  (II.B.18)

y=γ _(i) ln(I _(i))+γ_(j) ln(I _(j))+γ₀ ln(I ₀)+γ_(d)  (II.B.19)

Equations (II.B.18) and (II.B.19) are two linear equations in ln(I_(i)),ln(I_(j)), and ln(I₀). The coefficients, usingΔ=2(x_(i)y_(j)−x_(j)y_(i)), are given by

$\begin{matrix}{{{\eta_{i} = \frac{y_{j}}{\beta\Delta}},{\eta_{j} = {- \frac{y_{i}}{\beta\Delta}}},{\eta_{0} = {- \frac{\left( {y_{j} - y_{i}} \right)}{\beta\Delta}}},{\eta_{d} = \frac{\left( {y_{j} - y_{i}} \right)d^{2}}{\Delta}}}{{\gamma_{i} = {- \frac{x_{j}}{\beta\Delta}}},{\gamma_{j} = \frac{x_{i}}{\beta\Delta}},{\gamma_{0} = \frac{\left( {x_{j} - x_{i}} \right)}{\beta\Delta}},{\gamma_{d} = {- \frac{\left( {y_{j} - y_{i}} \right)d^{2}}{\Delta}}}}} & \left( {{{II}.B}{.20}} \right)\end{matrix}$

For the hexagonally packed sensor array in FIG. 1, with j≡i mod 6, andi=1, 2, . . . ,6, one has:

$\begin{matrix}{{{x_{i} = {d\mspace{14mu} {\cos \left( {i\frac{\pi}{3}} \right)}}},{x_{j} = {d\mspace{14mu} {\cos \left\lbrack {\left( {i + 1} \right)\frac{\pi}{3}} \right\rbrack}}}}{{y_{i} = {d\mspace{14mu} {\sin \left( {i\frac{\pi}{3}} \right)}}},{y_{j} = {d\mspace{14mu} {\sin \left\lbrack {\left( {i + 1} \right)\frac{\pi}{3}} \right\rbrack}}}}} & \left( {{{II}.B}{.21}} \right)\end{matrix}$

Substituting these into (II.B.20) gives

${\eta_{i} = \frac{\sin \left\lbrack {\left( {i + 1} \right)\frac{\pi}{3}} \right\rbrack}{\beta \; d\sqrt{3}}},{\eta_{j} = {- \frac{\sin \left( {i\frac{\pi}{3}} \right)}{\beta \; d\sqrt{3}}}},{\eta_{0} = {- \frac{\cos \left\lbrack {\left( {{2i} + 1} \right)\frac{\pi}{6}} \right\rbrack}{\beta \; d\sqrt{3}}}},{\eta_{d} = \frac{{\cos \left\lbrack {\left( {{2i} + 1} \right)\frac{\pi}{6}} \right\rbrack}d}{\sqrt{3}}}$${\gamma_{i} = {- \frac{\cos \left\lbrack {\left( {i + 1} \right)\frac{\pi}{3}} \right\rbrack}{\beta \; d\sqrt{3}}}},{\gamma_{j} = \frac{\cos \left( {i\frac{\pi}{3}} \right)}{\beta \; d\sqrt{3}}},{\gamma_{0} = \frac{\sin \left\lbrack {\left( {{2i} + 1} \right)\frac{\pi}{6}} \right\rbrack}{\beta \; d\sqrt{3}}},{\gamma_{d} = {- \frac{{\sin \left\lbrack {\left( {{2i} + 1} \right)\frac{\pi}{6}} \right\rbrack}d^{2}}{\sqrt{3}}}}$

The coefficient values for the seven sensor hexagonal ommatidium of FIG.2B are given in Table 1 and Table 2 of in FIG. 5B and each have thedimension of length.

FIG. 5A shows a functional block diagram of an analog circuitimplementation of equations (II.B.18) and (II.B.19) that can be realizedby summing the logarithmic amplification of the sensor signals with thecoefficients of Table 1 and Table 2 of in FIG. 5B. A first advantage ofthis embodiment is that the accuracy of the calculations depends only onthe accuracy of the logarithmic amplifier and the other analogcircuitry, and does not have discretization inaccuracies. Thecalculations also avoid inaccuracies due to time-interval averagingeffects that occur when using an analog-to-digital converter to samplethe sensor signals at fixed sample times. The sensors are denoted byR_(k) for k from 0 to 6, with R₀ denoting the central sensor. Theconfiguration 510 shows that the output of the sensor R₀ is received bya logarithmic amplifier, denoted ln, to produce an amplified outputsignal S₀. This same configuration is also shown being applied to thesensor outputs of sensors R_(k) for k from 1 to 6. The sigma symbols, asconventional, indicate summing operations, and the coefficientsdisclosed above are stored as constant sources.

The circuit yields six pairs of solutions 512 as given by the equationsabove.

While implementation using analog circuitry is a preferred embodiment,the digital sampling of the sensor responses could still be used toimplement the solutions given in equations (II.B.18) and (II.B.19). Foran ommatidium with an alternative configuration, a comparable analysiswill yield an alternate set of coefficients.

That the method for this configuration yields six pairs of solutionsimplies that up to and including six independently moving objects,producing six separate image features on the ommatidium can beseparately located. When the differencing detection method is appliedover real time, the separate detections allow for tracking the sixmoving objects. Generalizing from this to other configurations ofommatidia further implies that if there are a total of n many sensors inan ommatidium with n−1 arranged around a central sensor, then such anommatidium could be able to detect and track up to n−1 objects.

C. Orthogonal Decomposition Methods and Systems

The methods and systems described in this section can provide location,motion tracking, velocity and acceleration of an image featuretraversing the ommatidium. As with the differencing methods and systems,in preferred embodiments the processing may use real time parallelanalog processing. Other embodiments may process the sensor signalsusing sampling with digital circuitry. Exemplary methods generate anorthogonal decomposition of the ommatidium sensor output responses togenerate continuous position signals and, through temporal derivativesof the position signals, continuous velocity and acceleration signals ofthe image feature traversing the ommatidium. As with previous methodsand systems, hyperacuity can be achieved.

There are two aspects of the embodiments directed to the orthogonaldecomposition methods. The first aspect is related to the creation ofthe orthogonal signals. The second aspect relates to the processing ofthe orthogonal signals to generate a radial distance (rho: φ andincident angle (phi: ϕ, sometimes also denoted φ) with respect to theommatidium center of an image feature being tracked as it traverses theommatidium. At both stages, velocity and acceleration of the imagefeature can be generated by taking the temporal derivatives of thelocation/position signals.

As an overview, in some embodiments output signal intensity of theindividual ommatidium sensors are correlated with an orthogonal systemto provide location and motion tracking of an image feature traversingthe ommatidium field of view. The sensor outputs of the ommatidiumsensors are weighted with sine and cosine functions and summed using analgorithm to affix position to the sensors within the ommatidium toprovide a position of an image feature with respect to the center of theommatidium. The sums of the weighted sensor outputs represent the realand imaginary mathematical parts of the orthogonal decomposition of theimage feature position. For a moving objects, orthogonal intensitycomponents are constantly changing. The ratio of these values can beused to calculate the direction of object motion. Therefore, a change intrajectory can be extracted in real time. The temporal derivative of theposition signal yields the velocity of the image feature. The secondderivative yields the image feature's acceleration.

Embodiments of the orthogonal decomposition methods and systems can beunderstood with reference to the particular embodiment shown in FIG. 6Afor the seven sensor hexagonal ommatidium configuration of FIG. 2B. Asdescribed above, the sensors are arranged geometrically to form anommatidium disposed in a hexagonal pattern. Each sensor has a unique andwell-defined location relative to the central sensor which is also thecentral or reference point of the ommatidium. A beam of light incidenton this array of sensors will form an image-spot on a given sensor andon the extension of the Gaussian profiles of neighboring sensors. Thoughthe following description is with respect to that configuration, one ofskill in the art will see how the methods and systems can be applied toommatidia of other configurations.

FIG. 6A shows the (x, y) coordinate system for the ommatidium as awhole. In this configuration the angle formed from by the rays from theorigin to the center of two adjoining sensors is π/3. FIG. 6A also showsthe transformation matrix that can be used to convert the (x, y)coordinate system to a rotated coordinate system (x′, y′) for a givenangle of rotation γ. The method can be applied in the rotated coordinatesystem.

The orthogonal decomposition of the ommatidium's sensor signals can bebest described by means of complex functions and reference to FIG. 6B.In the case here, special optical coupling in each sensor results in aGaussian sensor response profile which is spatially dependent on theradial distance measured from the center axis of the sensor. ThisGaussian profile can be expressed, for the kth-sensor, as

$\begin{matrix}{{{s_{k}\left( {{\hat{I}}_{k},r_{k}} \right)} = {{\eta (\lambda)}{\hat{I}}_{k}e^{{- \beta}\; r_{k}^{2}}}},{k = 0},1,\ldots \mspace{14mu},6} & \left( {{{II}.C}{.1}} \right)\end{matrix}$

In (II.C.1), s_(k)(Î_(k),r) is the sensor response to an assumed thinbeam of light of intensity Î_(k) illuminating the kth-photoreceptorsystem normally at a distance r_(k) from the center of the sensor, η(λ)is an efficiency factor that includes the photon efficiency of thesensor, the optical transmissivity of the lens system, and other losseswhich are wavelength (λ) dependent. The signals from the sensor ontowhich the light beam is impinging and that of the neighboring sensorscan be added up to give a resultant system output signal S shown in FIG.6B. A simulated response of the sensors is as shown in FIG. 4 and wasdescribed above.

Significantly more useful information can be obtained by weighting eachsignal s_(k) by the complex exponential

$\begin{matrix}{W_{k} = {e^{{ik}\frac{\pi}{3}} = {{{\cos \left( {k\frac{\pi}{3}} \right)} + {i\mspace{14mu} {\sin \left( {k\frac{\pi}{3}} \right)}}} = {w_{k} + {i\mspace{14mu} {\overset{\_}{w}}_{k}}}}}} & \left( {{{II}.C}{.2}} \right)\end{matrix}$

This complex function relates the angular positions of the sensors tothe (x, y) frame of reference of the sensors arrangement or assembly forone ommatidium. These weights are indicated in FIG. 6A. Aftermultiplying (II.C.1) and (II.C.2) and summing over all k one obtains thecomposite complex signal

$\begin{matrix}{{U + {iV}} = {\frac{\sum\limits_{k = 1}^{6}{{\hat{I}}_{k}e^{{- \beta}\; r_{k}^{2}}e^{{ik}\frac{\pi}{3}}}}{\sum\limits_{k = 0}^{6}{{\hat{I}}_{k}e^{{- \beta}\; r_{k}^{2}}}} = \left\{ {\frac{\sum\limits_{k = 1}^{6}{{\hat{I}}_{k}e^{{- \beta}\; r_{k}^{2}}{\cos \left( {k\frac{\pi}{3}} \right)}}}{\sum\limits_{k = 0}^{6}{{\hat{I}}_{k}e^{{- \beta}\; r_{k}^{2}}}} + {i\frac{\sum\limits_{k = 1}^{6}{{\hat{I}}_{k}e^{{- \beta}\; r_{k}^{2}}{\sin \left( {k\frac{\pi}{3}} \right)}}}{\sum\limits_{k = 0}^{6}{{\hat{I}}_{k}e^{{- \beta}\; r_{k}^{2}}}}}} \right\}}} & \left( {{{II}.C}{.3}} \right)\end{matrix}$

It is convenient to call U the real part or component, and V theimaginary part or component of this decomposition of the sensor signals,both of which are mutually orthogonal to each other. FIG. 6B shows acircuit component diagram that can be used to calculate U and V. Asabove, the sensors are denoted R₀ . . . R₆ as shown producing signals s₀. . . s₆. The weights are stored as part of the cosine multiplier blocksw₀ . . . w₆ and the further sine blocks. The weights for the alternate(x′, y′) coordinates can be modified based on an input for the givenangle of rotation γ.

Based on (II.C.3), one may calculate the magnitude and the phase of thecomplex signal sum and express the composite complex signal as

$\begin{matrix}{{U + {iV}} = {{\left( {U^{2} + V^{2}} \right)^{1/2} \cdot \exp}\left\{ {i\mspace{14mu} {\tan^{- 1}\left( \frac{V}{U} \right)}} \right\}}} & \left( {{{II}.C}{.4}} \right)\end{matrix}$

That is, one obtains the magnitude (U²+V²)^(1/2) and phase φ=tan⁻¹(V/U)which can be used to locate the position of the image light spot in theommatidium.

It should be noted that the U and the V signals described here comprisethe orthogonal components of the intensity of the ommatidium sensorsignals and should not be confused with the I, Q phase decompositionthat is utilized in communications and antenna systems. Even though themathematics of the orthogonal decomposition is comparable, theunderlying principle of using sensor intensity instead of phase isnovel.

For locating the light spot relative to the ommatidium image space, boththe angle φ subtended between the line joining the spot and theommatidium center and the x-axis of the ommatidium reference frame, andthe distance ρ from the center of the ommatidium to the spot locationmust be computed. The distance is computed as the two-norm of thenormalized-by-Σ orthogonal signals U and V, scaled by a number k_(σ)that depends on the one-sigma (or beta-value) of the Gaussian curveutilized. That is,

$\begin{matrix}\begin{matrix}{\rho = {k_{\sigma} \cdot r \cdot \left\{ {(U)^{2} + (V)^{2}} \right\}^{1/2}}} \\{= {k_{\sigma} \cdot r \cdot \left\{ {\left\lbrack \frac{\sum\limits_{k = 1}^{6}{{\hat{I}}_{k}e^{{- \beta}\; r_{k}^{2}}{\cos \left( {k\frac{\pi}{3}} \right)}}}{\sum\limits_{k = 0}^{6}{{\hat{I}}_{k}e^{{- \beta}\; r_{k}^{2}}}} \right\rbrack^{2} +} \right.}} \\\left. \left\lbrack \frac{\sum\limits_{k = 1}^{6}{{\hat{I}}_{k}e^{{- \beta}\; r_{k}^{2}}{\sin \left( {k\frac{\pi}{3}} \right)}}}{\sum\limits_{k = 0}^{6}{{\hat{I}}_{k}e^{{- \beta}\; r_{k}^{2}}}} \right\rbrack^{2} \right\}^{1/2}\end{matrix} & \left( {{{II}.C}{.5}} \right)\end{matrix}$

The coefficient r in equation (II.C.5) can be taken as the one-sigmavalue of the Gaussian response profile of each sensor.

For larger area of an image light spot, the sensor signal is theresponse corresponding to the integral, over the spot's cross-sectionalarea, of the Gaussian profile. This is described next with the help ofFIG. 7A. In FIG. 7A an image feature light spot of radius r₀ is atlocation (x₀, y₀) of the ommatidium, over a sensor with a Gaussianresponse profile 710. The center line of the image light spot is at adistance d from the ith-sensor's center located at (x_(i), y_(i))relative to the ommatidium reference frame.

The sensor response can now be computed as the integral of the Gaussianprofile over the region R over the sensor surface 720 that the lightspot covers. Thus, for the ith-sensor,

$\begin{matrix}{{S_{i} = {\underset{R}{\int\int}{\eta (\lambda)}{\hat{I}}_{i}e^{{- \beta}\; r_{k}^{2}}{ds}}},{i = 0},1,\ldots \mspace{14mu},6} & \left( {{{II}.C}{.6}} \right)\end{matrix}$

The integration region R consists of the circle defined by the imagelight spot of radius r₀. The integration of the sensor's Gaussianresponse profile is obtained by Riemann sums. That is, the sum of theproducts of the sensor Gaussian value g(ρ_(k)) at distance ρ_(k) withrespect to the ith-sensor center (x_(i), y_(i)), times the annularregion areas ΔS_(k) over the light-spot region R. A three-dimensionaldepiction of the annular slices used for the numerical integration isdepicted in the following FIG. 7C.

The integration is carried over by summing up the individual products ofthe sensor response (Gaussian profile) at distance ρ_(k) from thesensor's center and the surface area ΔS_(k) of the annular strip ofwidth Δρ_(k) bounded by the light spot circle. The distance between thecenter (x_(i), y₁) of the sensor and the center (x₀,y₀) of the imagelight spot is given by d. Letting the arc length of the central fiber ofthe annular region be a_(k), the area of the annular strip within theregion R of the light-spot is given by

ΔS _(k) =a _(k)·Δρ=ρ_(k)·(2θ_(kc))·Δρ  (II.C.7)

In (II.C.7), θ_(kc) is the angle subtended by the radial vector ρ_(k)when the tip of this vector touches either point P or P′ in FIGS. 3 and4, where P and P′ are the intersections of the central annular sectorfiber (dotted line in the annular region in FIGS. 7B and 7C) with thecircle boundary. It can be shown from the geometry of the problem, usingthe law of cosines, that this angle is given by

$\begin{matrix}{\theta_{kc} = {\cos^{- 1}\left\{ \frac{\rho_{k}^{2} + d^{2} - r_{0}^{2}}{2\rho_{k}d} \right\}}} & \left( {{{II}.C}{.8}} \right)\end{matrix}$

The product of the Gaussian value g(ρ_(k)) and the annular region areaΔS_(k) is the “volume” V_(k) given by

V _(k) =g(ρ_(k))·Δs _(k)  (II.C.9)

Substituting (II.C.7) into (II.C.9) yields

V _(k) =g(ρ_(k))·ρ_(k)·(2θ_(kc))·Δρ.  (II.C.10)

Summing all these and replacing (II.C.8) into (II.C.10) gives

$\begin{matrix}{V = {2 \cdot {\sum\limits_{k = 0}^{n}{{{g\left( \rho_{k} \right)} \cdot \cos^{- 1}}{\left\{ \frac{\rho_{k}^{2} + d^{2} - r_{0}^{2}}{2\rho_{k}d} \right\} \cdot \rho_{k} \cdot {\Delta\rho}}}}}} & \left( {{{II}.C}{.11}} \right)\end{matrix}$

The sum extends to all annular strips that fit in the region R fromρ_(k)=ρ₀=d−r₀ to ρ_(k)=ρ_(n)=d+r₀. The number n of annular stripsbetween these limits is controlled by the width Δρ of the annular stripsthat is selected for accuracy considerations. That is, the number ofstrips and the strip width are related as

$\begin{matrix}{{{\Delta\rho} = \frac{2 \cdot r_{0}}{n}},{n \in Z^{+}}} & \left( {{{II}.C}{.12}} \right)\end{matrix}$

Substituting (II.C.12) and the normalized Gaussian response function

$\begin{matrix}{{g\; \left( \rho_{k} \right)} = {{\frac{1}{2{\pi\sigma}^{2}} \cdot \exp}\left\{ {- \frac{\rho_{k}^{2}}{2\sigma^{2}}} \right\}}} & \left( {{{II}.C}{.13}} \right)\end{matrix}$

into (II.C.11) yields

$\begin{matrix}{V = {2 \cdot {\sum\limits_{k = 0}^{n}{{\frac{1}{2{\pi\sigma}^{2}} \cdot \exp}{\left\{ {- \frac{\rho_{k}^{2}}{2\sigma^{2}}} \right\} \cdot \cos^{- 1}}{\left\{ \frac{\rho_{k}^{2} + d^{2} - r_{0}^{2}}{2\rho_{k}d} \right\} \cdot \rho_{k} \cdot \frac{2 \cdot r_{0}}{n}}}}}} & \left( {{{II}.C}{.14}} \right)\end{matrix}$

Substituting

$\rho_{k} = {d - r_{0} + {\left( {k + 1} \right) \cdot {\Delta\rho}} + \frac{\Delta\rho}{2}}$

into (II.C.14) and performing simplifications gives

$\begin{matrix}{V = {\frac{2 \cdot r_{0}}{n \cdot \pi \cdot \sigma^{2}} \cdot {\sum\limits_{k = 0}^{n}{\exp {\left\{ {- \frac{\rho_{k}^{2}}{2\sigma^{2}}} \right\} \cdot \left( {d - r_{0} + {\left( {k - 1} \right) \cdot {\Delta\rho}} + \frac{\Delta\rho}{2}} \right) \cdot \cos^{- 1}} \left\{ \frac{\rho_{k}^{2} + d^{2} - r_{0}^{2}}{2\rho_{k}d} \right\}}}}} & \left( {{{II}.C}{.15}} \right)\end{matrix}$

This Riemann sum can be computed numerically in simulation software.

Referring now to FIG. 8, the top trace shows the overlap volume(mathematical convolution) of a circle shaped target (as in FIG. 7A)moving across the Gaussian Response profile of a photo sensor and theresulting signal output of that circle moving across the equator of thatportion of the focal plane of the ommatidia. The bottom trace shows thepercentage overlap of that same circular shaped target with adjacentphoto sensors whose Gaussian response overlaps with the sensor where thetarget is transiting an equatorial line.

The orthogonal decomposition method can be adapted to ommatidia withother configurations.

D. Tracking and Trajectory Determinations

This section describes systems and methods for continuous and analogtracking of the position, velocity and acceleration of an image featureusing real time parallel analog processing of orthogonally decomposedommatidium sensor outputs into continuous signals representing theradial distance and incident angle relative to the ommatidium referenceaxes of an image feature traversing the ommatidium's field of view andthe temporal derivatives of these position signals to provide continuousvelocity and acceleration signals of the image feature.

Further analog processing of the orthogonal signals U and V enables thedetermination and continuous tracking of the radial distance (ρ) and theradial angle (φ) of an image feature relative to the ommatidiumreference which is the center of the central sensor (R₀) in theommatidium using analog circuitry (see FIG. 8). Temporal derivatives ofρ and φ signals can provide both velocity and acceleration of the objecton a real time continuous basis.

The methods and systems will be described for the particular embodimentof the seven sensor hexagonal ommatidium configuration, but it will beclear to one skill in the art that the methods and systems can beadapted for other configurations. The values of U, V, ρ and φ and theirtemporal derivatives may also be obtained by digitizing the ommatidium'sanalog sensor outputs and numerically computing these values so theembodiments also cover this approach too, as would be apparent to one ofskill in the art. The derivatives may also be obtained by sampling theommatidium's output at discrete time instants.

FIG. 9 shows the seven sensor hexagonal ommatidium configuration 900, onwhich an image feature 910 is moving in time across various sensors ofommatidia along the trajectory 920. The image feature is presumedinitially to be smaller than an individual sensor 221. The U and Voutputs determine the radial distance (ρ) and the radial angle (φ) inreal time.

FIG. 10 depicts a functional block diagram 1000 for an ommatidium O 1010comprising seven sensors R_(i), as shown. The sensor outputs s_(i), i=0,1, . . . , 6 are weighted in the Weighting Block 1020 to produce thereal (U) and imaginary (V) orthogonal components and the sum S of allthe sensor signals. U and V signals are subsequently processed by analogcircuitry in the Functions Block 1030 to produce ρ and φ of the imagefeature spot relative to the center of the ommatidium in polarcoordinates, and also produce the corresponding first and secondderivatives to yield velocity and acceleration of the image featurespot.

This orthogonal decomposition enables the determination of the angle φand radial distance ρ of an image feature traversing the ommatidiumwhere ρ is, ρ=k_(σ)·r·√{square root over (U²+V²)} where r is the sensorradius. The angle φ is determined by U and V, by the formula

$\phi = {\tan^{- 1}{\frac{V}{U}.}}$

Thus, with real-time analog processing, from the two continuous signalsU and V the image feature position within the ommatidium is specified.

FIG. 11 shows a component level block diagram for computing the firstand second derivative values of ρ and φ from the U and V inputs. Thespeed of the computation is only limited by the speed of the circuitcomponents, and not by a sampling speed.

Referring now to FIG. 12A and FIG. 12B, a trajectory is shown of animage feature smaller than an individual detector traversing theommatidium. FIG. 12A illustrates the geometry used for determining atrajectory of the image feature by first methods based on twoobservations of φ and ρ at different time instants. FIG. 12B illustratesthe geometry used for determining a trajectory of the image feature bysecond method based on continuous observations of φ and ρ.

There are two methods indicated in FIG. 12A for discrete sampledeterminations of the derivatives of φ and ρ. First, shown as (a) inFIG. 12A, is based on two observations, (φ₁,ρ₁) and (φ₂,ρ₂), made by thesystem at two different times, t₁ and t₂, respectively, as the imagefeature traverses the array. Denoting the points observed by P₁=P(t₁)and P₂=P(t₂), the trajectory displacement S and angle η relative to thearray reference line or x-axis are determined from

S=[ρ₁ ²+ρ₂ ²−2ρ₁ρ₂ cos(φ₂−φ₁)]^(1/2)  (II.D.1)

and

$\begin{matrix}{\eta = {\phi_{2} - {\arcsin {\left\{ {\frac{\rho_{1}}{S}{\sin \left( {\phi_{2} - \phi_{1}} \right)}} \right\}.}}}} & \left( {{{II}.D}{.2}} \right)\end{matrix}$

Equation (II.D.1) is obtained by straightforward application of the lawof cosines, while (II.D.2) is obtained by applying the law of sines asfollows:

$\begin{matrix}{{\frac{\sin \left( {\phi_{2} - \phi_{2}} \right)}{S} = {\frac{\sin (\theta)}{\rho_{2}} = \frac{\sin (\alpha)}{\rho_{1}}}},} & \left( {{{II}.D}{.3}} \right)\end{matrix}$

from which

$\begin{matrix}{\alpha = {\arcsin {\left\{ {\frac{\rho_{1}}{S}{\sin \left( {\phi_{2} - \phi_{1}} \right)}} \right\}.}}} & \left( {{{II}.D}{.4}} \right)\end{matrix}$

From FIG. 12A (b), the trajectory angle η is obtained as η=φ₂−α, whichis equation (II.D.4).

An alternative method to determine the trajectory from the two (φ, φdata values is as follows: if

α=ρ₂ cos(φ₂)−ρ₁ cos(φ₁), b=ρ ₂ sin(φ₂)−ρ₁ sin(φ₁)  (II.D.5)

then

$\begin{matrix}{{S = \left( {a^{2} + b^{2}} \right)^{1/2}},{\eta = {\arctan \left( \frac{b}{a} \right)}}} & \left( {{{II}.D}{.6}} \right)\end{matrix}$

This method appears to be simpler to implement as it does not requireinverse trigonometric functions as in the Discrete Method-1.Implementation of these two methods can be achieved using analogcircuitry employing of a time-delay line, or digitally from two (φ, φsamples.

Referring to FIG. 12B for a continuous signal approach, deriving thetrajectory is based on calculus. Starting with the exaggerated anglesubtended between the two points P₁=P(t₁) and P₂=P(t₂), in FIG. 12B, onecan make the following limiting observations. First, the arc s isobtained as s=ρΔφ. Relating to FIG. 12A for the discrete case, here onehas ρ₁=ρ and ρ₂=ρ+Δρ. Second, as t₂−t₁→dt, the following happens: Δρ→dρso that ρ₂→ρ₁=ρ, s→ds≈ds′ (that is, the arc approaches the secant),Δφ→dφ, the angle ∠P₁QP₂→π/2, and S→dS. Therefore, for an infinitesimaltime interval dt, one has

dS≈(ds ² +dρ ²)^(1/2)=[(ρdφ)² +dρ ²]^(1/2)  (II.D.7)

from which

$\begin{matrix}{\frac{dS}{dt} = {{\overset{.}{S} \approx \frac{{\rho^{2}d\; \phi \frac{d\; \phi}{dt}} + {d\; \rho \frac{d\; \rho}{dt}}}{\left\lbrack {{\rho^{2}d\; \phi^{2}} + {d\; \rho^{2}}} \right\rbrack^{1/2}}} = {\frac{{\rho \; d\; \phi \frac{d\; \phi}{dt}} + {\frac{d\; \rho}{\rho}\frac{d\; \rho}{dt}}}{\left\lbrack {{d\; \phi^{2}} + \left( \frac{d\; \rho}{\rho} \right)^{2}} \right\rbrack^{1/2}}.}}} & \left( {{{II}.D}{.8}} \right)\end{matrix}$

Since for most cases (dρ/φ²≈0, then (II.D.8) can be written as

$\begin{matrix}{\frac{dS}{dt} = {{\overset{.}{S} \approx \frac{{\rho^{2}d\; \phi \frac{d\; \phi}{dt}} + {d\; \rho \frac{d\; \rho}{dt}}}{d\; \phi}} = {{\rho \frac{d\; \phi}{dt}} + {\frac{1}{\rho}\frac{d\; \rho}{d\; \phi}{\frac{d\; \rho}{dt}.}}}}} & \left( {{{II}.D}{.9}} \right)\end{matrix}$

One can write

${\frac{d\; \rho}{d\; \phi} = \frac{d\; {\rho/{dt}}}{d\; {\phi/{dt}}}},$

so (II.D.9) becomes

$\begin{matrix}{\frac{dS}{dt} = {\overset{.}{S} \approx {{\rho \frac{d\; \phi}{dt}} + {\frac{1}{\rho \left( {d\; {\phi/d}\; t} \right)}{\left( \frac{d\; \rho}{dt} \right)^{2}.}}}}} & \left( {{{II}.D}{.10}} \right)\end{matrix}$

It can be shown that all the time-derivatives in (II.D.10) can beobtained from the time derivatives of U and V, that is, from dU/dt anddV/dt, which are given by

$\begin{matrix}{\phi = {\left. {\arctan \left( \frac{V}{U} \right)}\Rightarrow\frac{d\; \phi}{dt} \right. = {\frac{{U\frac{dV}{dt}} - {V\frac{dU}{dt}}}{U^{2} + V^{2}}.}}} & \left( {{{II}.D}{.11}} \right) \\{\rho = {\left. {k_{\sigma}{r\left( {U^{2} + V^{2}} \right)}^{1/2}}\Rightarrow\frac{d\; \rho}{dt} \right. = {k_{\sigma}r{\frac{{U\frac{dU}{dt}} + {V\frac{dV}{dt}}}{\left( {U^{2} + V^{2}} \right)^{1/2}}.}}}} & \left( {{{II}.D}{.12}} \right)\end{matrix}$

Substituting (II.D.11) and (II.D.12) into (II.D.10) gives

$\begin{matrix}{\frac{dS}{dt} \approx {{\rho \cdot \frac{{U\frac{dV}{dt}} - {V\frac{dU}{dt}}}{U^{2} + V^{2}}} + {\frac{1}{\rho} \cdot \frac{1}{\frac{{U\frac{dV}{dt}} - {V\frac{dU}{dt}}}{U^{2} + V^{2}}} \cdot {\left( {k_{\sigma}r\frac{{U\frac{dU}{dt}} + {V\frac{dV}{dt}}}{\left( {U^{2} + V^{2}} \right)^{1/2}}} \right)^{2}.}}}} & \left( {{{II}.D}{.13}} \right)\end{matrix}$

Using ρ=k_(σ)r(U²+V²)^(1/2) and after some algebra one obtains

$\begin{matrix}{\frac{dS}{dt} \approx {\rho \cdot \left\{ \frac{\left( \frac{dU}{dt} \right)^{2} + \left( \frac{dV}{dt} \right)^{2}}{{U\frac{dV}{dt}} - {V\frac{dU}{dt}}} \right\}}} & \left( {{{II}.D}{.14}} \right)\end{matrix}$

This is the velocity of the image spot over the ommatidium image space.Integration of (II.D.14) over time would give displacement, and thederivative of (II.D.14) with respect to time would give theacceleration. These can be achieved using an integrator and adifferentiator components in an analog circuit.

The angle η(t) of the instantaneous trajectory shown in FIG. 12B (b) isobtained as follows:

$\begin{matrix}{\alpha = {{\arctan \left( \frac{d\; \rho}{ds} \right)} = {{\arctan \left( \frac{\frac{d\; \rho}{dt}}{\frac{ds}{dt}} \right)} = {\arctan \left( \frac{\frac{d\; \rho}{dt}}{\rho \frac{d\; \phi}{dt}} \right)}}}} & \left( {{{II}.D}{.15}} \right)\end{matrix}$

Substituting (II.D.11) and (II.D.12) into (II.D.15) and using

$\eta = {\alpha + \phi - \frac{\pi}{2}}$

yields

$\begin{matrix}{{\eta (t)} = {{\arctan \left( \frac{{U\frac{dU}{dt}} + {V\frac{dV}{dt}}}{{U\frac{dU}{dt}} - {V\frac{dV}{dt}}} \right)} + \phi - \frac{\pi}{2}}} & \left( {{{II}.D}{.16}} \right)\end{matrix}$

FIG. 13 shows a block diagram of a circuit for obtaining the outputs ofequations (II.D.14) and (II.D.16). When so implemented, the trajectoryof the image feature can be calculated in real time.

E. Dynamic Contrast Control

This section describes embodiments for systems and methods for real-timetracking of an image feature with respect to the background andautomatically adjusting to contrast reversals between an image featureand background. Such systems and methods can be used, for example, whenan ommatidium-based system is imaging or tracking a dark object across abright sky. Light sensors of the ommatidium can then have the bright skybackground imaged onto most of their surfaces, with only a relativelysmall dark image feature on the sensor. As some previous embodimentspresume a bright object image and dark background, a potential problemor inaccuracy can be prevented by applying an inversion-biasing to theoutput signals of the sensors.

The theory of operation for dynamic contrast reversal can be understoodwith reference to FIG. 14 and FIG. 15. In the case of orthogonaldecomposition methods and systems, tracking objects regardless ofwhether lighter or darker than the ambient background, may be based ondifferentiating the U, V, and/or S-signals and by checking the sign, (+)or (−), of the differentiators' response signals. The sign can thentrigger a command to another section of the circuit that applies theinversion-biasing mechanism, or to revert to the normal operation (i.e.,a bright object against a dark background). While the image featurechanges brightness relative to the background as it traverses the array,the differentiators and sign detectors will enable or disable theinversion-biasing algorithm automatically so that the tracking continuesusing the appropriate tracking mode.

FIG. 14 shows combinations that are possible between the object imagefeature and the background:

Differentiators Signal Levels & Image-Background Combinations OutputSigns Bright Image (BI), Dark Background (BI) > (DB) => (+) (DB) . . .Bright Image (BI), Bright Background (BI) > (BB) => (+) or (BI) < (BB) .. . (BB) => (−) Dark Image (DI), Dark background (DI) > (DB) => (+) or(DI) < (DB) . . . (DB) => (−) Dark Image (DI), Bright Background (DI) <(BB) => (−) (BB) . . .The first case is shown in 1410, in which the image feature 1412 isbrighter than the background on the sensor. The graph for 1410 shows across section of the sensor response I. The sensor outputs due to thebackground I_(Bk) and the “target image” I_(Ti) satisfy I_(Bk)<I_(Ti).Thereafter, in the third case 1420 in which the image feature 1422 isdarker than the background I_(Bk)>I_(Ti). The detection scheme willautomatically switch between one mode (normal) or the other mode(inversion and biasing) according to the relative intensities betweenthe object image feature and the background, as the object image featuretraverses the ommatidium field of view. When the situation of 1420 isdetected, the detection scheme applies inversion 1430 so that the imagefeature 1432 is brighter than the background, and follows this withbiasing 1440 so that both outputs are now positive and the image feature1442 is now brighter than the background.

FIG. 15 shows an embodiment of a functional block diagram of a circuitthat can implement the dynamic contrast reversal. The particularembodiment shown applies to the case of U/V/S signals obtained using anorthogonal decomposition method to the seven sensor hexagonal ommatidiumstructure shown in FIG. 2B. The differentiators activate a switch toenable or disable the following inversion-and-biasing components forcontinual tracking of the object image feature.

F. Combined Methods and Systems

This section describes methods and systems that combine previouslydescribed methods and systems of locating or tracking an image featureof an object in a single ommatidium. Combining the differencing methodsdescribed in section II.B with the orthogonal decomposition methods ofsection II.C can produce improved accuracy in position, velocity andacceleration determinations.

There are several analytical advantages gained by expressing theorthogonal and differencing methods in compact form using matrix-vectornotation besides clarity and notational economy. In what follows,

(η_(k)) refers to the orthogonal system using Lagrange's planarhexagonal packing index and number η_(h), and

(T₂) refers to the differencing method for the triangular number T₂structure embedded in η_(h) and T_(n) for n>2, i.e., T₂⊂η_(h). Thisallows the implementation of zooming algorithms Z_(n) ^(n+1) for bothsystems, as will be described below in section III.C. Another advantageis the possibility of generation of a multiplicity ofmathematically-derived systems from these two basic ones by inversion,complementation, by dual spaces, by application of bilineartransformations, etc.

The orthogonal system is described here by using complex exponentialnotation. For the differencing method, T₂ is defined to be the distancebetween sensor centers for any T_(n) system besides the structure itrefers to.

For the orthogonal decomposition system, denoted

(η_(k)), one has

$\begin{matrix}{{{{{\overset{\_}{\vartheta}}_{\bot}\left( \eta_{h} \right)} = {\begin{bmatrix}\rho \\\phi\end{bmatrix} = \begin{bmatrix}{{{\overset{\_}{S}}^{T}\overset{\sim}{E}}} \\{\angle \left( {{\overset{\_}{S}}^{T}\overset{\sim}{E}} \right)}\end{bmatrix}}},{where}}{{\overset{\_}{S} = {{\overset{\_}{S}(t)} = \left\lbrack {s_{i}(t)} \right\rbrack}},{\overset{\sim}{E} = \left\lbrack {\exp \left\{ {j\; \overset{\_}{\Theta}} \right\}} \right\rbrack},{\overset{\_}{\Theta} = \left\lbrack {\theta_{i}\left( \eta_{h} \right)} \right\rbrack},{U = {\Re \; e\left\{ {{\overset{\_}{S}}^{T}\overset{\sim}{E}} \right\}}},{V = {\; m\left\{ {{\overset{\_}{S}}^{T}\overset{\sim}{E}} \right\}}},{{{and}\mspace{14mu}\sum} = {\overset{\_}{S}}_{1}}}} & \left( {{{II}.F}{.1}} \right)\end{matrix}$

For the differencing system, denoted

(T₂), one has

$\begin{matrix}{{{\overset{\_}{\vartheta_{\Delta}}\left( T_{2} \right)} = \left\lbrack p_{k} \right\rbrack},} & \left( {{{II}.F}{.2}} \right) \\{{{\overset{\_}{p}}_{k} = {\begin{bmatrix}x \\y\end{bmatrix}_{k} = {\begin{bmatrix}{\overset{\_}{\eta}}^{T} \\{\overset{\_}{\gamma}}^{T}\end{bmatrix}_{k}\begin{bmatrix}\overset{\_}{L} \\1\end{bmatrix}}_{k}}},} & \left( {{{III}.F}{.3}} \right)\end{matrix}$

where

${\begin{bmatrix}\overset{\_}{L} \\1\end{bmatrix}_{k} = \left\lbrack {{\ln \left( s_{0} \right)},{\ln \left( s_{k} \right)},{\ln \left( s_{k + 1} \right)},1} \right\rbrack_{k}^{T}},{\begin{bmatrix}\eta^{T} \\\gamma^{T}\end{bmatrix}_{k} = \begin{bmatrix}{\eta_{0,k},{\eta_{k,}\eta_{k + 1}},\eta_{d,k}} \\{\gamma_{0,k},\gamma_{k},\gamma_{k + 1},\gamma_{d,k}}\end{bmatrix}},{and}$${\eta_{0,k} = {- \frac{y_{k + 1} - y_{k}}{\beta\Delta}}},{\eta_{k} = \frac{y_{k + 1}}{\beta\Delta}},{\eta_{k + 1} = {- \frac{y_{k}}{\beta\Delta}}},{{\eta_{d,k} = \frac{\left( {y_{k + 1} - y_{k}} \right)T_{2}^{2}}{\Delta}};{y_{i} = {T_{2}{\sin \left( {\frac{\pi}{3}i} \right)}}}}$${\gamma_{0,k} = \frac{x_{k + 1} - x_{k}}{\beta\Delta}},{\gamma_{k} = {- \frac{x_{k + 1}}{\beta\Delta}}},{\gamma_{k + 1} = \frac{x_{k}}{\beta\Delta}},{{\gamma_{d,k} = {- \frac{\left( {x_{k + 1} - x_{k}} \right)T_{2}^{2}}{\Delta}}};{x_{i} = {T_{2}{\cos \left( {\frac{\pi}{3}i} \right)}}}}$for β = 4ln (2)/T₂², Δ = 2(x_(k)y_(k + 1) − x_(k + 1)y_(k)).

The methods may be combined using proportional-integral-derivative (PID)methods:

$\begin{matrix}{{\overset{\_}{\xi} = {\begin{bmatrix}x_{opt} \\y_{opt}\end{bmatrix} = {\underset{``P"}{f\left( {\rho,\phi,{\overset{\_}{p}}_{k}} \right)} + {A \cdot \underset{``I"}{\left\{ {\overset{R}{\underset{:P}{\Uparrow}}{{\overset{\_}{\vartheta}}_{\bot}\left( \eta_{h} \right)}} \right\}}} + \underset{``D"}{\left( {1 - A} \right) \cdot \left\{ {\Gamma {{\overset{\_}{\vartheta}}_{\Delta}\left( T_{2} \right)}} \right\}}}}},} & \left( {{{II}.F}{.4}} \right)\end{matrix}$

where f(ρ,φ,p _(k)) is a function of its arguments, being the P part inthe PID combining approach, which may include any correction terms,pre-acquisition coordinates as in the case of inter-ommatidia statepropagation, etc., and where

${\overset{R}{\underset{:P}{\Uparrow}}\begin{bmatrix}\rho \\\phi\end{bmatrix}} = {{\rho \cdot \begin{bmatrix}{\cos (\phi)} \\{\sin (\phi)}\end{bmatrix}} = \begin{bmatrix}x \\y\end{bmatrix}}$

represents the polar to rectangular conversion operator, and

${\Gamma = {\delta {\left\{ {\left( {{{\overset{R}{\underset{:P}{\Uparrow}}{{\overset{\_}{\vartheta}}_{\bot}\left( \eta_{h} \right)}}} - {{\overset{\_}{p}}_{k}}} \right) \cdot w} \right\} \cdot I_{6}}}},$

wherein δ is the Dirac delta function, and Γ is a diagonal kth-vectorSelect Matrix.

The matrix A=w·I₂ is a 2 by 2 diagonal weight matrix, where 0≤w≤1 wouldbe determined on specific tracking criteria and for particular cases. Ifw is set to 1, then the Orthogonal Decomposition Method would be usedexclusively. On the other extreme, if w is set to 0, then theDifferencing Method would be used exclusively. Intermediate values of wcan be set by different tracking strategies or criteria. For example, ifa large object is being observed, then the weight w could be set by theratio of the object size or angular extent divided by the ommatidia FOVextent. Other criteria may involve contrasts, brightness, object speed,distance from the ommatidia tracking limit or from object's distancefrom the ommatidia center (where the differencing method will have allthe positions converging). One can also set the weight w based onprobability measures obtained from the tracking data, from a feedbackloop of the tracking platform, or from an optimization trackingalgorithm that compares the two methods alternatively by switching theweight w between 0 and 1 and compares the object's trajectory and noise,etc. In summary, the algorithm(s) for determining of the optimal valuefor the combining weight w would be defined for each specificapplication. If w is set to zero, then also Γ=I₆ and multiple objectscan be tracked simultaneously.

These combined methods and systems can be implemented in analog ordigital circuitry, and may provide improved location detection andtracking of one or more independent image features.

III. Multiple Ommatidia Systems and Methods

This section describes methods and systems that use one or more arraysof sensors that have been organized into subarrays, with each subarrayacting as a single ommatidium, to perform detection and/or tracking ofone or more objects.

A. Multiple Ommatidia Systems

The methods and systems disclosed for a single ommatidia can be appliedwhen the sensor array comprises more than one individual ommatidia. Themethods for single ommatidia can be extended and adapted to combine thedetection and tracking results of each constituent ommatidium.Characteristics of the object or of multiple objects can be tracked,including size, shape, brightness, relative contrast which is directlyrelated to object detectability, and the rate of the object(s)traversing the field of view of the system.

Various embodiments of configurations that combine multiple ommatidia,each comprising multiple sensors, are shown in FIG. 16A, FIG. 16B, andFIG. 16.C, and described in further detail below. These figures showthree embodiments: FIG. 16A is a planar arrangement of seven ommatidiaunder a single focusing lens; FIG. 16. B comprises three ommatidia withindependent focusing lenses (one lens per ommatidium); and FIG. 16Cshows a particular example of numerous ommatidia arranged over ahemispherical (now a three-dimensional) supporting structure. It will berecognized that the configuration of FIG. 16C is a specific embodimentof embodiments in which the supporting structure is convex hull, or anarbitrary three-dimensional body in general. Each individual, or unit,ommatidium of a multiple ommatidia system may have any one of theconfigurations shown in FIGS. 2A to 2D, or another configuration ofmultiple sensors configured to implement hyperacuity.

Methods the can be used in multiple ommatidia systems, such as todetermine image feature position, include: orthogonal decompositionmethods, differencing methods, or the combination methods describedpreviously. Other methods will be clear to one of skill in the art. Thecombination orthogonal-differencing methods can use weighted methods,such as by a PID.

FIG. 16A shows a first configuration 1610, which may be implemented asan integrated circuit. Each individual ommatidium has the configurationof the seven sensor hexagonal ommatidium previously described. Theseindividual ommatidium are then themselves configured as six ommatidiahexagonally arranged around a central ommatidium. The counter clockwisenumbering of the sensors within each individual ommatidium is given bythe second index number; the first index indicates the ommatidiumnumber. The identification of the sensors within each ommatidiaarrangement will always be given by ordered pairs (o, s), where odenotes ommatidia number in an array, including a single ommatidium, ands denotes the sensor number within the oth-ommatidia. On each array, thecenter of the sensor (0, 0) would be taken as the origin of the array,which will be considered as the origin of the global reference frame.Likewise, the center of the sensor (o, 0) for the oth-ommatidia would beconsidered as the origin of the local reference frame for the particularommatidia.

A similar alternate configuration 1620 uses circular shaped sensors. Thecentral sensor of each seven sensor hexagonal ommatidium is theindicated with a bold circle.

In FIG. 16B, the embodiment configuration 1630 comprises multipleindependent ommatidia 1631 and 1632, each containing its own focusinglens, and put together mechanically on a planar support. In this casethe sensors are of different ommatidia are not adjacent, and larger gapsexist between the three shown in this arrangement. The arrangement 1630can be altered to have the three ommatidia touch each other as threetangent circles 120 degrees apart, as in configuration 1640. An “arraycenter” can be defined for both such configurations.

The configuration 1640 of FIG. 16B can function as a hierarchicalversion of the configuration shown in FIG. 2C.

In FIG. 16C, the embodiment configuration 1650 is a 3-D hemisphericalarrangement of unit ommatidia. As the other arrangements, a global arraycenter and frame relative to which all the ommatidia observations arereferred can also be defined as the center of the sphere the hemisphere.This is described in detail below.

FIG. 17A shows the geometry 1710 used for the methods for theconfiguration 1620 of FIG. 16A. The centers (x_(i), y_(i)) are given inTable 3 shown in FIG. 17C. Note that the radial distance of the circlescontaining all the six peripheral unit ommatidia centers is d√{squareroot over (f)}7, for d the radius of each circular sensor. The numbersare the cosines (for the X-coordinates) and the sines (for theY-coordinates) of the respective angles θ_(i) shown in Column 1 of Table3. These angles are the ones subtended from the reference line (thepositive part of the X-axis) and the reference point (the center of thecentral ommatidium) to each peripheral unit ommatidia in the arraysshown in FIG. 17A.

FIG. 17B shows the transformation 1720 or mapping of the position vectorfrom a local unit ommatidium coordinate frame to the global coordinateframe is obtained by the following formulas in polar coordinates:

$\begin{matrix}{\rho_{k} = \left\{ {\left\lbrack {x_{oi} + {\rho_{i}{\cos \left( \phi_{i} \right)}}} \right\rbrack^{2} + \left\lbrack {y_{oi} + {\rho_{i}{\sin \left( \phi_{i} \right)}}} \right\rbrack^{2}} \right\}^{1/2}} & \left( {{{III}.A}{.1}} \right) \\{\phi_{k} = {\arctan \left\{ \frac{y_{oi} + {\rho_{i}{\sin \left( \phi_{i} \right)}}}{x_{oi} + {\rho_{i}{\cos \left( \phi_{i} \right)}}} \right\}}} & \left( {{{III}.A}{.2}} \right)\end{matrix}$

A method for the transformation from a local to the global frame can beimplemented using the following pseudo-code as follows:

(A) Read the (x_(i),y_(i)) values from Table 3, or compute using:

-   -   (A.1) θ₀₁=tan⁻¹(√{square root over (3)}/2)    -   (A.2) For i=2 to 6

{ (A.3) θ_(0i) = θ₀₁ + (i − 1) · 60 (A.4) x_(0i) = d · √{square rootover (7)} · cos(θ_(0i)) (A.5) y_(0i) = d · √{square root over (7)} ·sin(θ_(0i))  }

-   -   (B) Read the ommatidium ρ_(i), φ_(i) values for the actual        spot-image location    -   (C) a=x_(0i)+ρ_(i) cos(φ_(i))    -   (D) b=y_(0i)+ρ_(i) sin(φ_(i))    -   (E) ρ_(k)=√{square root over (a²+b²)}    -   (F) θ_(k)=tan⁻¹(b/a) (Note: This can be accomplished with a        version of the arctan function that accounts for the quadrant of        (a, b).)

FIG. 18 shows the geometry 1800 for the configuration of multi-ommatidiacovering a body, such as shown in FIG. 16C. Once the local ommatidiavectors have been transformed to the global array vectors, thetransformation of the k^(th)-array global frame vectors onto a mainsystem reference frame of a convex hull or hemisphere containing all thearrays, can be subsequently applied as described next, using FIG. 18.The embodiment shown in FIG. 16C is for a hemispherical body for theommatidia arrays.

FIG. 18 shows a vector transformation from a local unit ommatidia arrayof sensors 1820 on a body 1810 to global coordinates for the body as awhole. The global coordinate system for the spherical coordinates shownin FIG. 18 is an orthogonal rectangular (x, y, z) frame whose origin andz-axis coincide with the center of symmetry and axis of symmetry of thehemisphere, respectively. The polar angle Θ_(i) for thei^(t)-photoreceptor or terminator is measured relative to the z-axis,and the azimuthal angle Φ_(i) is measured counter-clockwise from thex-axis of the right-handed reference frame. The x-axis is thepreferential or reference axis, such as for a supporting structure orvehicle. Each photo sensor (in the embodiment that a single ommatidiumis at each location) or each ommatidia array (e.g. 49-sensor arrays) inthe hemisphere has its own local reference plane and an associated localframe that can be defined as follows.

The normal-to-the-surface unit vector to the i^(th)-photoreceptor orarray is given, as can be derived from FIG. 18 relative to the system's(x, y, z) reference frame, as it is parallel to the radial vector R_(i),by:

û _(i) =Λu _(x) ,u _(y) ,u _(z)λ=Λ sin(Θ_(i))·cos(Φ_(i)),sin(Θ_(i))·sin(Φ_(i)), cos(Θ_(i))λ  (III.A.3)

This vector defines also the local horizon plane, a plane tangent to thehemispherical surface at the center of the ommatidium, or at theommatidia array's center. The south-looking unit vector ŝ_(i) isperpendicular to û_(i) and is contained in the plane defined by û_(i)and the z-axis (the trace or intersection between this plane and thehemisphere is called the local meridian of east-longitude Φ_(i)). Thislocal-frame unit-vector (or versor), ŝ_(i), is defined by

ŝ _(i) =

s _(x) ,s _(y) ,s _(z)

=

cos(Θ_(i))·cos(Φ_(i)), cos(Θ_(i))·sin(Φ_(i))m−sin(Θ_(i))

  (III.A.4)

Finally, the east-looking unit-vector ê_(i) forms with the other two aright-handed system. Hence, it is defined by

ê _(i) =û _(i) ×ŝ _(i) =

e _(x) ,e _(y) ,e _(z)

=

sin(Φ_(i)), cos(Φ_(i)),0

  (III.A.5)

The local frame given by the set {û_(i), ŝ_(i), ê_(i)} for thei^(th)-ommatidia, or array of ommatidia, is required for mapping anyobject's trajectory with respect to this local frame first, and thentransform or map the image coordinates (or image position vectors) tothe system or hemisphere frame given by the unit direction vectors{î_(i), ĵ_(i), {circumflex over (k)}_(i)}. This system frame is assumedto be the body (the convex hull hemisphere of FIG. 18) which supportsthis system of ommatidia. So the hemisphere's image position vectorswould in turn need to be transformed onto the body's reference frame bymeans of a chain product of homogeneous coordinates transformationmatrices. Such matrices may include scaling, rotations, and translationsfrom one frame to the next. An example of a homogeneous transformationmatrix is given by

$\begin{matrix}{\begin{pmatrix}x_{0} \\y_{0} \\z_{0} \\1\end{pmatrix} = {\begin{pmatrix}R_{11} & R_{12} & R_{13} & {\Delta \; x} \\R_{21} & R_{22} & R_{23} & {\Delta \; y} \\R_{31} & R_{32} & R_{33} & {\Delta \; z} \\0 & 0 & 0 & 1\end{pmatrix} \cdot \begin{pmatrix}x_{i} \\y_{i} \\z_{i} \\1\end{pmatrix}}} & \left( {{{III}.A}{.6}} \right)\end{matrix}$

In more compact form, equation (III.A.6) can be written as

$\begin{matrix}{\begin{pmatrix}{\overset{\rightarrow}{x}}_{0} \\1\end{pmatrix} = {\begin{pmatrix}R & {\Delta \; \overset{\_}{x}} \\0 & 1\end{pmatrix} \cdot \begin{pmatrix}{\overset{\rightarrow}{x}}_{i} \\1\end{pmatrix}}} & \left( {{{III}.A}{.7}} \right)\end{matrix}$

In (III.A.7), {right arrow over (x)}₀ is the position vector of theimage spot relative to the global frame, {right arrow over (x)}_(i) isthe position of the image spot relative to the ith-ommatidia, R is acomposite 3-axes Euler rotation matrix given by, R=Π_(k=1)³R_(k)(α_(k)), and Δx is a displacement vector representing thetranslation of the ommatidia or ommatidia array origin to the globalframe. The rotation matrix R can be written as

$\begin{matrix}{R = {{\prod\limits_{k = 1}^{3}{R_{k}\left( \alpha_{k} \right)}} = {{\begin{bmatrix}1 & 0 & 0 \\0 & {\cos \left( \alpha_{x} \right)} & {- {\sin \left( \alpha_{x} \right)}} \\0 & {\sin \left( \alpha_{x} \right)} & {\cos \left( \alpha_{x} \right)}\end{bmatrix}\begin{bmatrix}{\cos \left( \alpha_{y} \right)} & 0 & {\sin \left( \alpha_{y} \right)} \\0 & 1 & 0 \\{- {\sin \left( \alpha_{y} \right)}} & 0 & {\cos \left( \alpha_{y} \right)}\end{bmatrix}}{\quad\begin{bmatrix}{\cos \left( \alpha_{2} \right)} & {- {\sin \left( \alpha_{z} \right)}} & 0 \\{\sin \left( \alpha_{z} \right)} & {\cos \left( \alpha_{z} \right)} & 0 \\0 & 0 & 1\end{bmatrix}}}}} & \left( {{{III}.A}{.8}} \right)\end{matrix}$

Multiple ommatidia systems can provide improved detection and trackingcapabilities, as will now be described. Other embodiments andapplications include providing zooming detection and wider fields ofview, as will be disclosed in sections III.C and section IV.

B. Inter-Ommatidia Feature Tracking Methods and Systems

In detection and tracking using multiple ommatidia systems, it may occurthat an image feature either is near a boundary between two (or more)constituent unit ommatidia. Methods are now described for maintainingdetection and tracking with hyperacuity.

FIG. 19 shows the geometry 1900 for a method of tracking a moving imagefeature in a multi-ommatidia system as it transitions between theommatidia in the system's array. The particular embodiment is describedfor the 7 hexagonal planar ommatidia configured as described in FIG.16A, 1620. However, it can be understood that the methods to bedescribed can be adapted to other configurations.

During the transition, that is, when the light spot of the image featurebegins to appear on a neighboring ommatidia and begins to fade out ofthe present ommatidia, the local vectors would point to the center ofmass of that portion of the image light spot that affect the particularommatidia. Because of this, the transformed vectors of each ommatidiaexcited by the light spot would differ by some relatively small amount.If precision tracking is required, there are two methods that can befollowed to obtain the best estimates for the global set (ρ_(k),φ_(k)):(i) Perform a weighted average of the two sets obtained by thetransformations (1) and (2), where the weights would be computed as theratio of the intensities (equivalent to areas) of the parts of the spoton each ommatidia, divided the total light spot intensity. That is, if Iis the total light spot intensity, and I₁ and I₂ are the respective spotintensities detected on each ommatidia, then the weights w₁=I₁/I, andw₂=I₂/I, respectively. Then, these weights are used to perform aweighted average in both rho's and thetas obtained in the twotransformations; or (ii), utilize a dynamic state estimator/propagator(or state predictor).

FIG. 22 shows a flow chart for an exemplary method 2200 that adaptspreviously described detection and tracking methods for multipleommatidia systems. The steps are as disclosed previously but have nowadded the step 2210 which uses the size of the derivatives of U and V todetermine if motion is occurring. The test is whether the sum of theabsolute value of the derivatives exceeds a specified threshold ε. Instep 2220, in the case that the radial distance is found to exceed aminimum threshold, indicating the image feature is near anotherommatidium, a kinematic state predictor determines whether anotherommatidium has also detected the image feature. In the case, tested atthe Om On decision branch, that another ommatidium indeed has alsodetected the image feature the method proceeds to either compute anoverlap or continue with sensor readings.

C. Meta-Ommatidia Zooming Methods

This section describes methods and systems by which, in a multipleommatidia system, the separate ommatidia can be used hierarchically toprovide location zooming. In some of these embodiments a locationobtained for an image feature with respect to global coordinates of themultiple ommatidia system is used to select a component ommatidium ofthe ommatidia system, and the signals of that component ommatidium arethen used to provide a finer resolution location. These hierarchicalconfigurations of multiple ommatidia may be generalized to include theability to zoom using both the differencing and the orthogonaldecomposition methods previously described.

First embodiments of zooming make use of the orthogonal decompositionapproach of section II.C. FIG. 20 illustrates how a multiple ommatidiasystem 2010 with 49 sensors can have its signals summed to produce areduced set of seven signals and so effectively act as a seven sensorsingle ommatidium 2020. This can allow for less processing when theeffective ommatidium is providing desired detection and trackingaccuracy.

FIG. 21 shows an alternate configuration 2110 of a nine sensorommatidium. The particular arrangement into triangles of 2110 allows forsensor signal summation to produce a reduced signal set to function asthe reduced ommatidium 2120.

D. Dynamic Reconfiguration In Multiple Ommatidia Systems

In the multiple ommatidia systems and methods described so far, it waspresumed that the individual sensors of the entire system were uniquelyassigned to respective unit ommatidia. This section describes systemsand methods by which the arrangement of the individual sensors of theentire system can be dynamically reconfigured “on the fly” to provideimproved tracking and detection. Such methods and systems are describedwith respect to the particular embodiment shown in FIG. 23. However, itwill be clear to one of skill that the methods and systems can beadapted to other configurations.

FIG. 23 shows the configuration 2300 of sensors of the multipleommatidia configuration 1610, with internal sensors explicitlyenumerated. Also indicated are boundary sensors 2310. The methods maycomprise reselecting which individual sensors are to form the respectivesensors of unit seven sensor hexagonal ommatidia. Once, e.g., the sensorchosen to be the central sensor of the central ommatidium is selected,the reorganization of the sensors into respective unit ommatidia is thenchosen accordingly.

FIG. 24 shows a block diagram of the operations that perform thereorganization of the sensors. One limitation that may be enforced isthat only interior sensors can be selected as centers for unitommatidia. The methods and systems illustrated in FIG. 23 and FIG. 24may be used in the systems next described for multiple ommatidia systemsand methods for increased field of view.

IV. HYPRIS System Applications of Ommatidia Systems

The following describes further embodiments (referred to as HYPRIS™)that extend and apply the embodiments disclosed above. HYPRIS™embodiments provide continuous kinematic target tracking and may beexpanded to any arbitrary field of view by combining multiple ommatidiainto a multi-layered hierarchy, much like a compound eye structure.Inter-ommatidia processing and tracking may be accomplished using analogor sampling with numerical processing. In this section the objectsdiscussed above will also be referred to as targets.

HYPRIS™ embodiments are based on modularity of the multiple ommatidiamethods and systems. This modularity facilitates the placement ofommatidia on a smooth curved surface (such as vehicle hulls, aircraftwings, or missiles) providing full-surround situational awarenesswithout breaking aerodynamic contours. The ability to operate on awingtip rather than an airframe centerline may be made possible by realtime compensation for wing motion relative to the airframe. The compoundeye structure is useful for small or large airframes as it does notrequire pan and tilt control or large volume aplanatic optics.

A. Wide FoV Methods with Multiple Ommatidia

HYPRIS's modular system is extensible up to a full spherical FoV. Abasic single ommatidium module may operate over a limited field of view(FoV); in some case this may be approximately ten degrees. Parallel,analog processing of a single ommatidium produces the continuous analogsignals representing the kinematic states of up to six targets for theseven sensor ommatidium embodiment used as examples above. (Recall alsothat a greater number of targets can be tracked with otherconfigurations.) Two embodiments for creating a wide FoV using anommatidium are as follows. The first method assumes that the kinematicstate signals of each ommatidium, arranged on a compound surface toachieve the desired FoV, will be numerically sampled and processed usingalgorithms to facilitate inter-ommatidia tracking of targets across theentire FoV. Although numerical processing is employed to accomplishmulti-segment tracking, HYPRIS may use embedded processors to greatlyreduce the computational burden required by image-based kinematicanalysis of targets.

A second approach uses the multi-layer, hierarchical processingarchitecture possible with multi-ommatidia systems. This replicates thesame sensing and processing methods of a basic ommatidia module atsuccessively higher meta levels in order to preserve the ability todetect and track multiple targets over a wide FoV on a continuous analogbasis. Information is propagated from higher to lower layers to detect,isolate and track objects with increasingly greater spatial resolutionall in the analog domain. Since HYPRIS' s parallel analog processing isnot constrained by the Nyquist digital sampling rate and the processinglimitations of image based sensor systems, it allows for very high speedoperation.

The methods and systems may implement the capacity to calculate andsubtract self-motion of a imaging sensor platform already exists butrequires a separate reference system such as an INS or GPS for absoluteaccuracy or a large computational burden if self-motion is estimatedfrom image-based analysis by the sensor system. In further embodiments,the methods and systems may implement embodiments that can be used toestimate and extract self-motion using kinematic state signals, whichalready contain motion cues, with algorithms and processing at a greatlyreduced computational burden.

B. IR and Polarization Detection

In some embodiments, the multi-ommatidia configurations may use IR orpolarization sensitive sensors.

C. High Level Partition Methods

FIG. 25 shows a block diagram of a multiple ommatidia system asdescribed in relation to FIG. 16A. The multiple ommatidia system isshown, in this embodiment, using an orthogonal decomposition methoddescribed in section II.C.

FIG. 26 shows a block diagram 2600 of how the multiple ommatidia systemof FIG. 25 can be implemented to use both differencing and orthogonaldecomposition methods. The outputs of the differencing method applied tothe unit ommatidia comprise the locations, velocities and accelerations,2610, with according to respective rectangular coordinates systems ofthe unit ommatidia. Further outputs 2622 and 2624 are obtained usingorthogonal decomposition methods, and comprise respective radialdistances and incident angles from each of the unit ommatidia. Theoutputs may be combined to produce net location and trajectoryinformation.

FIG. 27 shows a component level block diagram of the multiple ommatidiasystem of FIG. 26. As shown, it can be implemented with continuousreal-time analog components to avoid problems due to sampling limits andlatency.

By implementing one or more multiple ommatidia systems, HYPRIS methodsand systems can be configured to detect and track multiple targets overa wide field of view.

What is claimed is:
 1. A method of calculating a position of a feature of an object at hyperacuity accuracy, the method comprising: receiving a plurality of output signals from a plurality of sensors, wherein: the plurality of sensors have overlapping spatial response profiles that are non-linear and axio-symmetric about respective centers of each of the plurality of sensors; and the plurality of sensors produce the plurality of output signals based on the overlapping spatial response profiles; generating first weighted signals by weighting each of the plurality of output signals with respective first sinusoidal weights; generating second weighted signals by weighting each of the plurality of output signals with respective second sinusoidal weights; generating a first orthogonal resultant signal, wherein the first orthogonal resultant signal comprises a sum of the first weighted signals; generating a second orthogonal resultant signal, wherein the second orthogonal resultant signal comprises a sum of the second weighted signals; and determining the position of the feature using the first orthogonal resultant signal and the second orthogonal resultant signal.
 2. The method of claim 1, wherein the overlapping spatial response profiles comprise approximately Gaussian profiles.
 3. The method of claim 1, wherein the method is performed with analog circuitry in real time.
 4. The method of claim 1, further comprising using the first orthogonal resultant signal and the second orthogonal resultant signal to determine a radial distance and an incident angle to the position of the feature relative an origin of the plurality of sensors.
 5. The method of claim 4, further comprising determining at least one of a velocity or an acceleration of the position of the feature using at least one of: time derivatives of the radial distance and the incident angle; or time derivatives of the first orthogonal resultant signal and the second orthogonal resultant signal.
 6. The method of claim 4, wherein: the plurality of sensors comprises optical sensors arranged around a central optical sensor; and the origin of the plurality of sensors is a center of the central optical sensor.
 7. The method of claim of claim 1, wherein: the plurality of sensors comprises optical sensors arranged on a focal plane of an objective lens; and the feature is a component of an image of the object focused onto the focal plane by the objective lens.
 8. The method of claim 1, wherein: the first sinusoidal weights comprise cosines of successive multiples of π/3 radians; and the second sinusoidal weights comprise sines of successive multiples of π/3 radians.
 9. The method of claim 1, further comprising: determining a plurality of positions of the feature over time; and determining a trajectory of the feature using a plurality of positions of the feature.
 10. The method of claim 1, wherein the plurality of sensors are arranged in a hexagonal array comprising a center sensor and six sensors surrounding the center sensor.
 11. A system for determining a position of a feature of an object at hyperacuity accuracy comprising: a plurality of sensors that have overlapping spatial response profiles that are non-linear and axio-symmetric about respective centers of each of the plurality of sensors, wherein the plurality of sensors produce a plurality of output signals based on the overlapping spatial response profiles; processing circuitry configured to: generate first weighted signals by weighting each of the plurality of output signals with respective first sinusoidal weights; generate second weighted signals by weighting each of the plurality of output signals with respective second sinusoidal weights; generate a first orthogonal resultant signal, wherein the first orthogonal resultant signal comprises a sum of the first weighted signals; generate a second orthogonal resultant signal, wherein the second orthogonal resultant signal comprises a sum of the second weighted signals; and determine the position of the feature using the first orthogonal resultant signal and the second orthogonal resultant signal.
 12. The system of claim 11, wherein the overlapping spatial response profiles comprise approximately Gaussian profiles.
 13. The system of claim 11, wherein the processing circuitry comprises analog circuitry that operates in real time.
 14. The system of claim 11, wherein the processing circuitry is further configured to use the first orthogonal resultant signal and the second orthogonal resultant signal to determine a radial distance and an incident angle to the position of the feature relative an origin of the plurality of sensors.
 15. The system of claim 14, wherein the processing circuitry is further configured to determine at least one of a velocity or an acceleration of the position of the feature using at least one of: time derivatives of the radial distance and the incident angle; or time derivatives of the first orthogonal resultant signal and the second orthogonal resultant signal.
 16. The system of claim 14, wherein: the plurality of sensors comprises optical sensors arranged around a central optical sensor; and the origin of the plurality of sensors is a center of the central optical sensor.
 17. The system of claim 11, wherein: the plurality of sensors comprises optical sensors arranged on a focal plane of an objective lens; and the feature is a component of an image of the object focused onto the focal plane by the objective lens.
 18. The system of claim 11, wherein: the first sinusoidal weights comprise cosines of successive multiples of π/3 radians; and the second sinusoidal weights comprise sines of successive multiples of π/3 radians.
 19. The system of claim 11, wherein the processing circuitry is further configured to: determine a plurality of positions of the feature over time; and determine a trajectory of the feature using a plurality of positions of the feature.
 20. The system of claim 11, wherein the plurality of sensors are arranged in a hexagonal array comprising a center sensor and six sensors surrounding the center sensor. 